# Dependency: sim.results in "sim_results.Rdata"

library(xlsx)

# Load simulation output

load("sim_results.Rdata")

true.effect <- sim.results[[1]]
sim.results.long <- sim.results[-1]

J <- 1000

# Plot distribution of ATTs

avg.beta <- true.effect[["ATT"]]

# Figure 1

par(mfrow=c(3,2)) 
#par(oma=c(6,1,0,0),mar=c(0,1,7,3)) 

plot(density(sim.results[["vATT.weight"]][,1]), xlim = c(0.5,2.5), main = "Propensity Score Weighting \n(without regression adjustment)", cex.main = 0.95, ylim = c(0,6))
lines(density(sim.results[["vATT.wweight"]][,1]), lty = 2)
abline(v = avg.beta)
legend("topright", c("unweighted", "weighted"), bty = "n", lty = c(1, 2))

plot(density(sim.results[["vATT.weight.reg"]][,1]), xlim = c(0.5,2.5), main = "Propensity Score Weighting \n(with regression adjustment)", cex.main = 0.95, ylim = c(0,6))
lines(density(sim.results[["vATT.wweight.reg"]][,1]), lty = 2)
abline(v = avg.beta)
legend("topright", c("unweighted", "weighted"), bty = "n", lty = c(1, 2))

plot(density(sim.results[["vATT.nn"]][,1]), xlim = c(0.5,2.5), main = "Nearest Neighbor Matching \n(without regression adjustment)", cex.main = 0.95, ylim = c(0,6))
lines(density(sim.results[["vATT.wnn"]][,1]), lty = 2)
abline(v = avg.beta)
legend("topright", c("unweighted", "weighted"), bty = "n", lty = c(1, 2))

plot(density(sim.results[["vATT.nn.reg"]][,1]), xlim = c(0.5,2.5), main = "Nearest Neighbor Matching \n(with regression adjustment)", cex.main = 0.95, ylim = c(0,6))
lines(density(sim.results[["vATT.wnn.reg"]][,1]), lty = 2)
abline(v = avg.beta)
legend("topright", c("unweighted", "weighted"), bty = "n", lty = c(1, 2))

plot(density(sim.results[["vATT.s"]][,1]), xlim = c(0.5,2.5), main = "Subclassification Matching \n(without regression adjustment)", cex.main = 0.95, ylim = c(0,6))
lines(density(sim.results[["vATT.ws"]][,1]), lty = 2)
abline(v = avg.beta)
legend("topright", c("unweighted", "weighted"), bty = "n", lty = c(1, 2))

plot(density(sim.results[["vATT.s.reg"]][,1]), xlim = c(0.5,2.5), main = "Subclassification Matching \n(with regression adjustment)", cex.main = 0.95, ylim = c(0,6))
lines(density(sim.results[["vATT.ws.reg"]][,1]), lty = 2)
abline(v = avg.beta)
legend("topright", c("unweighted", "weighted"), bty = "n", lty = c(1, 2))

### Calculate confidence interval at each j

for (i in 1:length(sim.results.long)) { # ITERATE OVER SIMULATIONS
  
  sim.results.long[[i]] <- cbind(sim.results.long[[i]], matrix(NA, nrow=J, ncol=3))
  
  sim.results.long[[i]][,3] <- sim.results.long[[i]][,1] - (1.96 * sim.results.long[[i]][,2])
  
  sim.results.long[[i]][,4] <- sim.results.long[[i]][,1] + (1.96 * sim.results.long[[i]][,2])
  
  sim.results.long[[i]][,5] <- ifelse(avg.beta <= sim.results.long[[i]][,4] & avg.beta >= sim.results.long[[i]][,3], 1, 0)
  
}

sim.results.array <- array(NA, c(length(sim.results.long), J, 5))
for (i in 1:length(sim.results.long)){
  sim.results.array[i,,] <- sim.results.long[[i]]
}

# Estimated treatment effects

est.mean <- apply(sim.results.array[,,1], 1, mean)
names(est.mean) <- names(sim.results.long)

# Bias

bias <- function(v,true) {
	sum(v - true)/length(v)
}

est.bias <- apply(sim.results.array[,,1], 1, bias, true = avg.beta)
names(est.bias) <- names(sim.results.long)

# MSE

mse <- function(v,true) {
	sum((v - true)^2)/length(v)
} 

est.mse <- apply(sim.results.array[,,1], 1, mse, true = avg.beta)
names(est.mse) <- names(sim.results.long)

# CI

est.ci <- t(apply(sim.results.array[,,1], 1, quantile, p = c(0.025,0.975)))
rownames(est.ci) <- names(sim.results.long)

# COVERAGE

est.cov <- apply(sim.results.array[,,5], 1, mean) * 100
names(est.cov) <- names(sim.results.long)

# Table 1

TABLE <- cbind(est.mean, est.bias, est.mse, est.ci, est.cov)

write.xlsx(TABLE, "sim_table.xlsx")

